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Triple-differential cross sections for two-photon double ionization of molecular hydrogen are pre- 
sented for a central photon energy of 30 eV. The calculations are based on a fully ab initio, 
nonperturbative, approach to the time-dependent Schrodinger equation in prolate spheroidal co- 
ordinates, discretized by a finite-element discrete-variable-representation. The wave function is 
propagated in time for a few femtoseconds using the short, iterative Lanczos method to study the 
correlated response of the two photoelectrons to short, intense laser radiation. The current results 
often lie in between those of Colgan et al [J. Phys. B 41 (2008) 121002] and Morales et al [J. Phys. 
B 41 (2009) 134013]. However, we argue that these individual predictions should not be compared 
directly to each other, but preferably to experimental data generated under well-defined conditions. 
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Mapping the four-body breakup processes of the hy- 
drogen molecule via multiphoton absorption has recently 
become feasible experimentally at the free-electron laser 
facility FLASH. Although the angular distributions and 
the kinetic-release-energy spectrum of ionic fragments of 
D2 were measured [H, determining the angular distribu- 
tion of the ejected electrons is very difficult. Compared 
to one-photon double ionization (DI) of H2 @, @| , an ac- 
curate theoretical description of the two-photon DI of 
the H 2 molecule is also extremely challenging. Previ- 
ous calculations by Colgan et al [J], who employed the 
time-dependent close-coupling (TDCC) method, and by 
Morales et al who used a time-independent exterior 
complex scaling (ECS) treatment, showed considerable 
disagreements regarding both the shape and the mag- 
nitude of the predicted triple differential cross section 
(TDCS). However, both predicted a strong dependence 
of the angular distribution of the ejected electrons on 
the alignment angle between the linear laser polarization 
vector (e) and the molecular axis (<;). 

The discrepancies between the previous, computation- 
ally very demanding calculations provided the primary 
motivation for the present project. We emphasize, how- 
ever, that these calculations were performed for much 
different conditions. Colgan et al [4] used a femtosecond 
laser pulse with a peak laser intensity of 10 15 W/cm 2 that 
could lead to serious depletion of the ground state and 
may not be in the perturbative regime. The calculation 
of Morales et al [f| was for a well-defined photon energy 
and is equivalent to a weak pulse of near infinite duration. 
In light of the ongoing discussions, it seemed important 
to perform another independent calculation to investi- 
gate the similarities and differences in the results and to 
understand their origin. Consequently, we also consider 
the two-photon DI of the H2 molecule in the parallel and 
perpendicular geometries at an incident photon energy 
centered around 30 eV and equal sharing of the excess 
energy. Our pulse (see details below) lies in the perturba- 
tive regime but, like the TDCC model, has a bandwidth 



that may access doubly excited electronic states in the 
H2 spectrum unavailable to Morales et al. The effect of 
such states on the outcome of a short-pulse experiment is 
not known at this time. The details could be highly com- 
plex, and a theoretical treatment may ultimately have to 
go beyond the current model. 

The two-center nature of H2 already destroys the 
spherical symmetry of the problem even in the absence 
of an external laser field. While it is numerically feasible 
to carry out the calculations in a spherical coordinate 
system at the center of the H2 molecule, this becomes 
more and more difficult for heavier diatomics due to the 
slow convergence of the multi-center electron-nuclear in- 
teraction. Fortunately, the prolate spheroidal coordinate 
system offers an attractive alternative for diatomics Q, 
since the electron-nuclear interaction is rendered benign. 

The addition of an external field causes a symmetry 
reduction when the laser field is not parallel to the inter- 
nuclear axis. Compared to the atomic He target, there- 
fore, the situation in this simplest two-electron molecule 
is complicated due to an additional degree of freedom, 
namely the alignment angle between the e- and <r-axes. 
While the emission modes of photoelectrons from atomic 
targets exhibit rotational symmetry with respect to the 
polarization vector, this may no longer be the case for an 
arbitrary e — <; geometry. All of these issues conspire to 
make the response of the H2 molecule to temporal laser 
fields much more computationally demanding than the 
two-photon DI of the helium atom. 

For the frozen-nuclei approximation used here and 
in [1, [|| to be physically meaningful, the time scale for 
ionization by the field has to be much shorter than the 
characteristic time of the nuclear motion. For direct, 
rather than sequential double ionization, this condition is 
fulfilled here, since the simultaneously ejected electrons 
have a speed over ten times larger than what the pro- 
tons can achieve in the subsequent Coulomb explosion. 
For a central photon energy of 30 eV, we expose the H2 
molecule to a 10-cycle sine-squared laser pulse at a peak 
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intensity of 10 14 W/cm 2 , followed by a 2-cycle field- free 
evolution. This "10 + 2"-cycle takes about 1.6 fs, which 
is sufficiently long to extract a well-defined cross section 
in the perturbative regime of relatively low peak intensi- 
ties. Nevertheless, a comparison with results from time- 
independent calculations is not directly possible, due to 
the energy width of several eV. 

Our treatment of the six-dimensional time-dependent 
Schrodinger equation (TDSE) of the laser-driven H 2 
molecule relies on the two-center prolate spheroidal co- 
ordinates (£,r),(p). The ranges of the variables are £ € 
[1, +oo), 77 e [—1, +1], and 99 6 [0, 2tt]. For the temporal 
laser field E(t) = eE(t), the TDSE in the dipole length 
gauge reads 
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2,t)= U1+U2 + — + E{t)-{r 1 +r 2 ) ¥(1,2, t). 
L r X2 i 

(1) 

Here Hi represents the field-free single-electron Hamil- 
tonian for electron "i" while is its coordinate, mea- 
sured relative to the center of the molecule. As usual, 
ru = |n - r 2 \. 

We expand the H2 wave function in the body-frame as 



777 1 7772 



n mim2 ^ lim ,^v2,t)e l{mivi+m2ip2 \ 

where mi and r«2 denote the magnetic quantum num- 
bers of the two electrons along the molecular axis. 
To discretize this partial differential equation, we em- 
ploy the finite-element discrete-variable-representation 
(FEDVR) approach for the (£,77) variables Q. Thus 
n mim2 (£i,?7i>£2!??2)i) is expanded in a product of "ra- 
dial" {/;(£)} and "angular" {g k (v)} DVR functions. The 
notation originates from the asymptotic behavior of £ — >• 
2r/R and rj — > cos 9 in spherical coordinates, when the 
electron is far away from the nuclei. 

To treat the electron-electron interaction we employ 
the Neumann expansion of I/V12 in prolate spheroidal 
coordinates |7j, 

Here £> and £< are the larger and smaller of £1 and £2, 
respectively. -Pj m '(£) and Q[ m '(£) are regular and irreg- 
ular Legendre functions defined for £ G (l,+oo), while 

Pj m '(r/) is specified for 77 6 [— 1,+1]. By suitably gen- 
eralizing the approach used in 8] for the spherical case, 
it is possible to reformulate the computation of the re- 
quired matrix elements of the operator in Eq. Q as the 
solution of a two-center Poisson equation. This results in 
a diagonal representation of the matrix elements of the 
electron-electron Coulomb interaction and considerably 
simplifies the FE-DVR discretization procedure. 

The time-dependent electron wave packets of the laser- 
driven H2 molecule were generated by using our recently 



developed Arnoldi-Lanczos algorithm on the DVR mesh 
points [9j . Compared to the treatment of atomic targets 
in spherical coordinates, however, the boundary condi- 
tions in spheroidal coordinates require further elabora- 
tion. As demonstrated in [6j , the "radial" and "angular" 
functions are not simple polynomials for odd \m\. Special 
attention must be paid to the square-root asymptotics as 
£ — > 1 (from above) and \r)\ — >■ 1 (from below) when |m| 
is odd. By introducing the factors -y/£ 2 — 1 and y/l — if , 
respectively, in the DVR bases /•»(£) and gk(v) f° r °dd 
|m| to account for the behavior near the boundary, it is 
possible to use a single quadrature for all m. 

A spatial box with £ max = 100 was set up to trun- 
cate the configuration space. Solving the TDSE in imag- 
inary time with a grid of iVg x N v = 224 x 9 and 
M max = |wi|max = |m 2 | m ax = 3, we obtained the elec- 
tronic energy of the initial X X E+ state at the internu- 
clear distance of R = 1.4 ao (ao = 0.529 x 10~ 10 m de- 
notes the Bohr radius) as —1.8873 atomic units (a.u.), 
in good agreement with the recent benchmark value of 
-1.888761428 a.u. by Sims and Hagstrom [l(J (after tak- 
ing out the nucleus- nucleus interaction of 1/1.4). This 
yields a double- ionization threshold of 51.4 eV above the 
initial electronic state. 

To calculate the angle-resolved differential cross sec- 
tions, we project the wave packet at the end of the time 
propagation onto uncorrelated two-electron continuum 
states, which are approximately constructed in terms of 
the single-electron continuum states of the H^~ ion. For a 
momentum k, the "radial" part of the continuum state, 
-^|raL(£)> behaves asymptotically as 
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Here Ai m i 9 (fe) is the two-center Coulomb phase shift, 
c = kR/2, a = 2R, and £ = \m\ + q for the Hj ion. In ad- 
dition to |m|, an integer q, which is the number of nodes 
of the angular function, is used to label the continuum 
state. In the unificd-atom limit, R — > 0, the above two- 
center radial function reduces to the well-known atomic 
Coulomb wave function. 

Let us comment again on the similarities and the dif- 
ferences between the various approaches. One might ex- 
pect our calculation to yield similar results to the TDCC 
model, although our numerical implementation is very 
different from that of Colgan et al [4|. For example, they 
used a finite-difference method and a flat-top pulse with 
12-cycle time duration at a much stronger intensity of 
10 15 W/cm 2 . Most importantly, however, both models 
project onto uncorrelated Coulomb functions. 

In the ECS approach of Morales et al ||, on the 
other hand, the information was extracted from a time- 
independent wave function, which effectively corresponds 
to an infinite propagation time. The complicated three- 
body correlated Coulomb boundary condition is avoided 
in this approach. We would best simulate the ECS model 
by exposing the H2 molecule to a weaker light field with a 
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longer pulse duration, and also by propagating for a long 
time after the external field has died off. Then the elec- 
trons have moved very far away from the nuclei and the 
projection to uncorrelated functions is increasingly ap- 
propriate [111 ]. Since this would require a large amount 
of computational resources, we decided on the compro- 
mise of testing our procedure by varying the number of 
field-free propagation cycles between one and three. We 
only noticed a small sensitivity to the time of projection, 
and hence are confident that our results are converged 
to a few percent. In addition, it should be noted that 
TDCC and ECS results for one-photon DI @, Q, per- 
formed with similar philosophies of extracting the infor- 
mation, showed excellent agreement with each other. 
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FIG. 1. (Color online) Energy probability distribution of two 
ejected electrons in the parallel (a) and perpendicular (b) ge- 
ometries for a sine-squared laser pulse of 10 optical cycles, a 
central energy of 30 eV, and a peak intensity of 10 14 W/cm 2 . 
The color bars correspond to multiples of 10~ 7 eV -2 . Note 
the different scales in the two panels. 

Figure [1] depicts the correlated energy distributions of 
the two ejected electrons after two-photon absorption. As 
they are similar in shape for the both parallel and the per- 
pendicular geometries, information about any possibly 
preferred direction of the outgoing electrons is smeared 
out in the energy distributions. Figure [T] reveals, how- 
ever, that the probability of two-photon DI in the per- 
pendicular geometry is about ten times larger than for 
the parallel case, in qualitative agreement with the find- 
ings of (H- 

Next we explore the angle-resolved cross sections in 
a particular plane, namely that formed by the molecular 
axis and the laser polarization axis, a coplanar configura- 
tion. Note that the angles of the two ejected electrons, 9\ 
and 6*2 , are measured with respect to the e vector rather 
than the <j axis. 

Our TDCS results for the parallel geometry are shown 
in Fig. [21 They did not change within the thickness of 
the lines when an enlarged spatial box of £ max = 150 
was used. They are numerically converged at |m| max = 5 
(see the convergence checks below). Near the maximum 
of the TDCS, the time- independent ECS results are sig- 
nificantly larger than the present FE-DVR and TDCC 
predictions. The present TDCS results for the paral- 
lel geometry are approximately three times smaller than 
those of Morales et al @ , but twice as large as those of 
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FIG. 2. (Color online) Coplanar TDCS for two-photon DI of 
H2 at equal energy sharing (Ei = E2 — 4.3 eV) of the two 
ejected electrons in the parallel geometry. The laser param- 
eters are the same as in Fig. [1] Also shown are the TDCC 
results of Colgan et al 0] and the ECS results of Morales et 
al [fl, multiplied by the scaling factors indicated in the leg- 
end. The definition of the TDCS is consistent in all three 
calculations. 



Colgan et al [4J, i.e., they generally lie in between the pre- 
vious predictions. This is indicated by the scaling factors 
shown in the legend. 

The various sets of results disagree in the predicted 
magnitude as well as in the shape of the angular de- 
pendence. A prominent difference concerns the "wing" 
structures of non vanishing TDCS values near 82 — 0° and 
360° for 6»i = 0° and 30° observed in Ref. 0]. In addi- 
tion to the dominant back-to-back escape mode, there is a 
noticeable forward emission in the prediction of Colgan 
et al. This counter-intuitive forward-to-forward escape 
mode for two electrons carrying the same kinetic energy 
and traveling in the same direction is not supported by 
either the ECS or the present FE-DVR calculations. 

Figure [3] displays the corresponding TDCSs for the 
perpendicular orientation. Here the situation is differ- 
ent from the previous case. Although the ECS results Q 
are generally larger than the present and the TDCC pre- 
dictions they are of similar magnitude. The agree- 
ment between the FE-DVR and TDCC calculations is 
quite satisfactory in this case. In contrast to the par- 
allel geometry, the shapes of the angular distributions 
in the perpendicular geometry are much closer to those 
obtained for the helium atom at 42 eV @. 

The TDCS for orientation of the molecular axis paral- 
lel to the polarization vector exhibits the molecular effect 
to the largest extent. Only the 1 S g — > 1 S U — > 1 S g path 
of transitions is open in this geometry, because the total 
magnetic quantum number along the <?-axis must be con- 
served. The final state, therefore, has the same symmetry 
as the initial state. At a right angle between the e- and 
axes, on the other hand, the two electrons may be ionized 
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FIG. 3. (Color online) Same as Fig. [2] for the perpendicular 
geometry. No scaling factors were applied to compare the 
various predictions. 
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FIG. 4. (Color online) Convergence of the TDCS predictions 
with increasing |m| max for the parallel (a) and perpendicu- 
lar (b) geometry. The results for |m| max = 4 and 5 are hard 
to distinguish. The laser parameters are the same as in Fig. [2] 



the selection rule here is \AM\ = 1 instead of \AM\ = 
(here M — vn\ + m,2 is the total magnetic quantum num- 
ber along the <r-axis). This is qualitatively similar to 
those of the atomic target through 1 S e —5- 1 P° — > 1 S e 
or 1 S e — > 1 P° — > 1 D e . An interference effect between 
the open channels 1 H g and 1 A g in H2 is evident, as it is 
in the helium atom between the 1 S e and 1 D e channels 
[l2| . Such interference cannot be observed in the parallel 
orientation, where only one channel is open. 



We also performed extensive convergence checks by in- 
creasing |m| max from 3 to 5. Figure [4] shows the con- 
vergence for the fixed electron angles of B\ = 0° and 
Q\ = 90°. In the parallel case, the TDCS calculated with 
|?7i|max — 3 at Q\ — 0° is already fully converged, but for 
Q\ = 90° |m| max — 5 is required to obtain essentially con- 
verged results. In the perpendicular case, on the hand, 
the TDCS at 0\ = 0° is very sensitive to the value of 



until it stabilizes at m 



5, while the TDCS 



at 6*i = 90° is less sensitive to the choice of |m| max . 

Although the parallel and perpendicular orientations 
exhibit different convergence patterns in terms of 9\ and 
02, we observe that the TDCSs in both geometries behave 
similarly when we transfer the result to the molecular 
frame. In that frame, when the fixed electron is measured 
in the direction perpendicular to the molecular axis, the 
back-to-back escape mode is most sensitive to the \m\ 
value. When one electron is detected along the molecular 
axis, the back-to-back mode is the most stable. 

In summary, we calculated the two-photon DI of the 
H2 molecule by an intense femtosecond laser pulse. Our 
TDCS results are generally closer to the predictions from 
the TDCC than from the ECS approach, except for the 
small-angle behavior in the parallel geometry, where we 
do not see what seems to be an unphysical increase of the 
TDCS. A magnitude difference of about a factor of two, 
rather than five [5(, remains in the parallel case. Once 
again, however, we emphasize that a straight compari- 
son of these results is inappropriate. Great care should 
be taken when such a comparison is made, particularly 
when experimental data become available. Knowing the 
details of the experimental setup will be essential. We 
plan to extend our present treatment to angular distri- 
butions in more general cases with an arbitrary alignment 
angle between the polarization vector and the molecular 
axis, and to further study the sensitivity of the theoreti- 
cal predictions on the details of the pulse. 
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